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Two particles, initially in a product state, become entangled when they come together and start to 
interact. Using semiclassical methods, we calculate the time evolution of the corresponding reduced 
density matrix pi, obtained by integrating out the degrees of freedom of one of the particles. To 
quantify the generation of entanglement, we calculate the purity V{t) — Tr[pi(t)^]. We find that en- 
tanglement generation sensitively depends (i) on the interaction potential, especially on its strength 
and range, and (ii) on the nature of the underlying classical dynamics. Under general statistical 
assumptions, and for short-ranged interaction potentials, we find that ■P(t) decays exponentially 
fast if the two particles are required to interact in a chaotic environment, whereas it decays only 
algebraically in a regular system. In the chaotic case, the decay rate is given by the golden rule 
spreading of one-particle states due to the two-particle coupling, but cannot exceed the system's 
Lyapunov exponent. 
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When two systems (. . . ) enter into temporary inter- 
action (■■■), and when after a time of mutual influence 
the systems separate again, then they can no longer he 
described in the same way as before, viz. by endowing 
each of them with a representative of its own. This is 
how entanglement was characterized by Schrodinger al- 
most seventy years ago [J] . Entanglement is arguably the 
most puzzling property of multipartite quantum systems, 
and often leads to counterintuitive predictions due to, in 
Einstein's words, spooky action at a distance. Entangle- 
ment has received a renewed, intense interest in recent 
years in the context of quantum information theory . 

In the spirit of Schrodinger's above formulation, one is 
naturally led to ask the following question: "What deter- 
mines the rate of entanglement production in a dynamical 
system ?" . Is this rate mostly determined by the interac- 
tion between two, initially unentangled particles, or does 
it depend on the underlying classical dynamics ? This is 
the question we address in this paper. Previous attempts 
to answer this question have mostly focused on numerical 
investigations 0, S IE IE S 111 > with claims that entan- 
glement is favored by classical chaos, both in the rate it 
is generated |E ^ Q and in the maximal amount it can 
reach jE|. In particular, strong numerical evidences have 
been given by Miller and Sarkar for an entanglement pro- 
duction rate given by the system's Lyapunov exponents 
Q . These findings have however been recently challenged 
by Tanaka et al. lEl, whose numerical findings show no 
increase of the entanglement production rate upon in- 
crease of the Lyapunov exponents in the strongly chaotic 
but weakly coupled regime, in agreement with their ana- 
lytical calculations relating the rate of entanglement pro- 
duction to classical time correlators. Ref. is seemingly 
in a paradoxical disagreement with the almost identical 
analytical approach of Ref. , where entanglement pro- 
duction was found to be faster in chaotic systems than in 
regular ones 9]. This controversy thus calls for a better 
analytical understanding of the problem. 

We present a semiclassical calculation of the time- 



evolved density matrix p{t) for two interacting, distin- 
guishable particles. In the usual way, entanglement is 
quantified by the properties of the reduced density ma- 
trix pi{i) = TT2[p(t)], obtained from the two-particle 
density matrix by tracing over the degrees of freedom 
of one (say, the second) particle. At t = 0, the two par- 
ticles are in a product state of two narrow wavepackets. 
We will quantify entanglement by calculating the purity 
■p(t) = Tr[pi(t)^] which varies from for fully entangled 
to 1 for factorizable two-particle states (the discus- 
sion of how this could translate into the violation of a Bell 
inequality llj is postponed to later works). Compared 
to other measures of entanglement such as the von Neu- 
mann entropy or the concurrence, the purity presents the 
advantage of being analytically tractable. For the weak 
coupling situation we are interested in here, numerical 
works have moreover shown that von Neumann and lin- 
ear entropy Sun = 1 — ■p(t) behave very similarly j^. 
We thus expect the purity to give a faithful and generic 
measure of entanglement. We note that our semiclas- 
sical approach is straightforwardly extended to the case 
of undistinguishable particles, provided the nonfactoriza- 
tion of the reduced density matrix due to particle statis- 
tics is properly taken care of 12] . 

Our approach is reminiscent of the semiclassical meth- 
ods developed by Jalabert and Pastawski 13, 14] (see 
also 1^^) in the context of the Loschmidt Echo. We will 
show how the off-diagonal matrix elements of pi are re- 
lated to classical action correlators, in a similar way as 
in Ref. [E ■ Under nonrestrictive statistical assump- 
tions, we find that, following an initial transient where 
pi relaxes but remains almost exactly pure, entangle- 
ment production is exponential in chaotic systems, while 
it is algebraic in regular systems. The asymptotic rate of 
entanglement production in chaotic systems depends on 
the strength of the interaction between the two particles, 
and is explicitely given by a classical time-correlator. We 
note that, as is the case for the Loschmidt Echo jl^lr^. 
this regime is also adequately captured by an approach 
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based on Random Matrix Theory (RMT) ^ - the time- 
correlator is then replaced by the golden rule spreading of 
one-particle states due to the interaction. For stronger 
coupling however, the dominant stationary phase solu- 
tion becomes interaction independent, and is determined 
only by the classical dynamics, the system's Lyapunov 
exponents giving an upper bound for the rate of entangle- 
ment production. The crossover between the two regimes 
occurs once the golden rule width becomes comparable to 
the system's Lyapunov exponent. Still, one has to keep 
in mind that long-ranged interaction potentials can lead 
to significant modifications of this picture, especially at 
short times, due to an anomalously slow vanishing of off- 
diagonal matrix elements of pi within a bandwidth set 
by the interaction range. 



We start with an initial two-particle product state 
iV'i) ® \ip2) = IV'ij'02)- The state of each particle is a 
Gaussian wavepacket ^/'i,2(y) = (ttct^)"^/** exp[ipo • (y — 
ri.2) — |y — ri,2p/2cr^]. We write the two-particle Hamil- 
tonian a.s TC = H^I + I(^H + U, i.e. the two particles 
are subjected to the same Hamiltonian dynamics. At 
this point, we only specify that the interaction potential 
U is smooth, depends only on the distance between the 
particles, and is characterized by a typical length scale ( 
(which can be its range, or the scale over which it fluc- 
tuates). Setting h = 1, the two-particle density matrix 
evolves according to p{t) — exp[— i7-^t]po exp[iHt] start- 
ing initially with po — V'2)(V'i; '/'2|- The elements 
/Oi(x,y;i) = J c?r(x, r|p(i)|y, r) of the reduced density 
matrix read 



pi(x,y;i) {na')-"! dv j ntidy, cxp[-{(yi - ri)^ + (y^ - ra)^ + (yg - ri)^ + (y4 - r2)2}/2a2] 

X exp[ipo(yi +y2 -ys -y4)] (x,r| exp[-im]lyi,y2)(y3,y4| exp[im]|y,r). (1) 
We next introduce the semiclassical two-particle propagator 

(x,r|exp[-im]|yi,y2) = {-iY ^Cl'J, eyi-p[i{Ss{yi,:>^]t) S'^' (y2 , r; -I- ^^^s' (yi, x; y2, r; - 

which is expressed as a sum over pairs of classical trajectories, labelled s and s', respectively connecting yi to x and 
y2 to r in the time t. Each such pair of paths gives a contribution containing one- (denoted by Ss and Ss') and two- 
particle (denoted by Ss^s' = Jq dtiU{qs{ti),<ls' {ti)) |2^) action integrals accumulated along s and s', a pair of Maslov 
indices ps and figi, and the determinant Cs^s' of the stability matrix corresponding to the two-particle dynamics in 
the (2(i)— dimensional space. (With the above definition, Cg.s' is real and positive.) We consider sufficiently smooth 
interaction potentials varying over a distance much larger than a. We thus set Sg.s' (yi, x; y2, r; t) ~ Ss^s' (ri, x; r2, r; t) 
(still we keep in mind that ri and r2, taken as arguments of the two-particle action integrals, have an uncertainty 
0(cr)), and use the narrowness of the initial wavepackets to linearize the one-particle actions in y^ — r j (i = 1, . . . 4; 
j = 1,2). We consider the weak coupling regime, where the one-particle actions vary faster than their two-particle 
counterpart. We thus perform a stationary phase approximation on the 5"s to get the semiclassical reduced density 
matrix as 

pi(x,y;i) = (~W)'^^ exp[i{5,(ri,x;t)~5,(ri,y;i)}] 

s,l 

X J dr^Ms.s'Ml exp[i{5s,s'(ri,x;r2,r;t) - 5i,s'(ri,y;r2,r;i)}], (2) 

s' 

Msy = Cl'l exp[-i-(^^ - exp[-y {(p^ - po)^ + (p^' - Po)^}]- 

It is straightforward to see that Tr[pi(i)] — 1 and /3i(x, y; t) — [pi(y, x; t)]*, as required. Enforcing a further stationary 
phase condition on Eq. (|2J) amounts to performing an average over different initial conditions ri 2. It results in s = s', 
x = y, and thus (/9i(x, y;t)) = (5x,y/^^ (^^ is the system's volume), i.e. only diagonal elements of the reduced density 
matrix have a nonvanishing average (the ergodicity of (/9i(x, x;t)} is due to the average over initial conditions). For 
each initial condition, p\ has however nonvanishing off-diagonal matrix elements, with a zero-centered distribution 
whose variance is given by (pi (x, y; t)/9i (y, x; i)). Squaring Eq. (jJJ, averaging over v\,2 and enforcing a stationary 
phase approximation on the S"s, one gets 

(pi(x,y;t)pi(y,x;t)) = {ina^f^J drrfr' ^ ^ Ms,s' M;,™ M^, MI„, (T), (3) 

T = exp[i{5s,s'(ri, x; r2, r; t) - Si^s' (ri, y; r2, r; t)}] exp[i{5i,m(ri, y; r2, r'; t) - 5s,m(ri, x; r2, r'; t)}]. (4) 



3 



Our analysis of Eqs. (|3I4|I starts by noting that is given by the sum of two positive contributions. First, those 

particular paths for which r = r' and s' = m, accumulate no phase = 1) and thus have to be considered separately. 
On average, their contribution does not depend on x nor y, and decays in time only because of their decreasing 
measure with respect to all the paths with r ^ r' . Their average contribution is given by 



dvidv2dr ^ \Ms,-m\^ \Mi 
.1 



oc 



^ exp[— (Ai + A2)i] ; chaotic, 
17-2 2 . regular. 



(5) 



This is obtained by taking half of the C's in the A^'s as Jacobians for the coordinate transformation from r's to p's, and 
replacing the other half by the large time limit Ts'l Cs,m = c'^^C^ oc (tQ/t)'^ for regular, and « exp[— (Ai + A2)i], for 
chaotic systems (In Eq. we explicitely wrote particle indices; in our case where the two particles are subjected to 
the same Hamiltonian, the Lyapunov exponents are equal, Ai = A2). Two facts are noticeable here: these contributions 
do not depend on the interaction strength, moreover, Eq. ((SJ gives a lower bound for the decay of (Ipil^). 

The second, generic contributions to (Ipil^) decay in time with (T). From Eq. it is natural to expect that (JF) 
is a decreasing function of |x — y| and t only. Sums and integrals in Eq. Q can then be performed separately to get 
(pi(x,y;i)pi(y,x;t)) = ^^(t) + (J^{^,y;t)) /Q^ , with 



(J^(x,y;i)) = exp 



2 / dhdt2 {U{(is{h),qs'(h))Uiqs{t2);(is'it2))) 



exp 



+ 2 / dt,dt2 {ll{(ls{tl),Cls'{h))U{cil{t2);(is'{t2))) 




(6) 



The behavior of (.F(x, y; t)) sensitively depends on the distance |x — y| between the endpoints of the classical trajec- 
tories q_s{t) and q;(t) connecting ri to x and y respectively. For small |x — y|, the second line on the right-hand side 
of Eq. almost compensates the first one (2^ , and a Taylor expansion of the difference of the two-particle action 
integrals in Eq.Q) gives in lowest order 



(^(|x-y| <C;0) = exp 



2 (x-y)4x-y)0 / dii di2(9i^Mq.(^i), Q^' (^i)) 4'^^(q«(^2); q.' (i2))) 



a, 13=1 



■(7) 



The reduced density matrix has a Gaussian decay away from the diagonal, which we expect to hold for |x — y| <C C 
(the distance over which the two terms in Eq. (|SJ almost cancel each other), and for short enough times t < t (after 
which very different, uncorrelated classical paths / and s may still reach two neighboring endpoints). An estimate for 
T is given by the time it takes for two initial conditions within a distance a (this is the uncertainty in ri) to move 
away a distance oc C from each other. In a chaotic system, this is the Ehrenfest time r — Aj~^ ln(C/(T), while in a 
regular system, r is much longer, r — 0{[(/a]^) 0. At larger distances, |x — y| :$> (, {T) is given by the first term 
on the right-hand side of Eq. Note that, because the four classical paths in that term come in two pairs, the 
dependence on |x — y| vanishes. 

This concludes our semiclassical calculation of the reduced density matrix for interacting two-particle systems. We 
have learned that the variance of off-diagonal matrix elements of pi is determined by classical correlators, with the 
important caveat that they are bound downward by the expressions given in Eq. Q. For the rest of the discussion, 
we note that, provided the correlators in Eqs. © and {Tj) decay faster than oc \ti — ^2!^^, the off-diagonal matrix 
elements exhibit a dominant exponential decay in time. This condition is rather nonrestrictive and is surely satisfied 
in a chaotic system. We therefore assume from now on a fast decay of the correlations. 



(W(q,(ti),q,.(ti))i^(q,(t2);q.'(i2))) = F (5(ti - ^2), 
di'^U{cis{ti),ci,,{h))d^fu{ci,{t2);cis-{t2))) = l5^,p5{ti-t2). 



(8) 
(9) 



Entanglement is quantified by the purity 'P(t) ~ 
J dxdy{pi{x,y;t)pi{y,x;t)), which is straightforward to 
compute from Eqs. (|3I9|I . We get three distinct regimes 
of decay: (i) an initial regime of classical relaxation for 
t < T. During that time, pi evolves from a pure, but 
localized pi(0) = |ri)(ri| to a less localized, but still al- 
most pure pi (t) , with an algebraic purity decay obtained 



from Eqs. and ® as V(t < r) ~ n-\7r/2-ftf/'^ 
(even in the case of a correlator saturating at a fi- 
nite value for \ti — ^2! ^ 00, which may occur in regular 
systems, this initial decay will still be algebraic oc t~'^); 
(ii) a regime where quantum coherence develops between 
the two particles so that pi becomes a mixture. From 
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Eq. lOl and the first line of Eq. © with Eq. © one gets 
T-'{t) ~ exp[— min(2r; Ai + X2)t] in a chaotic system; in a 
regular system one has V{t) oc {to/t)^^, since T,Q{t) dom- 
inates the decay of V{t) independently of the correlators 
in Eas. (|8l9l) : and (iii) a satm'ation regime where the pu- 
rity reaches its minimal value. In the chaotic limit, this 
saturation value can be estimated using a RMT approach 
as Vioo) = 2{a'^/n) + 0{n-^) ^J, which is in qualita- 
tive agreement with the results obtained in Ref.Q for 
the von Neumann entropy. There is no reason to expect 
a universal saturation value in the regular regime. 

Analyzing these results, we note that Eqs. © and 
(UI) are reminiscent of the results obtained for V{t) by 
perturbative treatments in Refs. 0, 01 , but they apply 
well beyond the linear response regime. Our weak cou- 
pling condition that the one-particle actions S vary faster 
than the two-particle actions S roughly gives an upper 
bound U < E for the typical interaction strength U {E 
is the one-particle energy). The linear response regime 
is however restricted by a much more stringent condition 
?7 < A <C (A is the mean level spacing) The 
decay regime (ii) of Pj^lreconciles the a priori contra- 
dicting claims of Refs. [MflH and Ref. j^. For weak cou- 
pling, the decay of V{t) is given by classical correlators, 
and thus depends on the interaction strength, in agree- 
ment with Ref. 0. However, V{t) cannot decay faster 
than the bound given in Eq. Q , so that at stronger cou- 
pling, and in the chaotic regime, one recovers the results 
of Ref. 3 . Simultaneously, regime (ii) also explains the 
data in Fig. 2 and 4 of Ref. (2i] , showing an exponential 



decay of 'P(t) in the chaotic regime, and a power-law de- 
cay with an exponent close to 2 in the regular regime (this 
power-law decay was left unexplained by the authors of 
Ref. 0)- Our semiclassical treatment thus presents a 
unified picture of the problem. 

Three more remarks are in order here. First the sta- 
tionary phase solutions leading to the above results still 
hold in the case when the two particles are subjected to 
different one-particle Hamiltonians. Second, the power- 
law decay of 'P{t) predicted above for regular systems, is 
to be taken as an average over initial conditions ri^2 (in 
that respect see Refs. Q and H), but may also hold 
for individual initial conditions, as e.g. in [3. Finally, 
there are cases when the correlators ((HJ and (PJ decay 
exponentially in time with a rate related to the spectrum 
of Lyapunov exponents. This also may induce a depen- 
dence of 'P{t) on the Lyapunov exponents, which can be 
captured by the linear response approach of Ref. . We 
note however that this is not a generic situation, as most 
fully chaotic, but nonhyperbolic systems have power-law 
decaying correlations. 

As a concluding line, and noting similarities between 
the problem treated here and that of decoherence by an 
environment |lO|. we anticipate that a semiclassical ap- 
proach as the one presented here (see also Ref. |2^) could 
clarify how decoherence relates to the intrinsic dynamics 
of the system ,26j . 
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